home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / ada / gwuada_9.zip / ARITH.C < prev    next >
C/C++ Source or Header  |  1993-07-27  |  47KB  |  2,081 lines
  1. /*
  2.  * Copyright (C) 1985-1992  New York University
  3.  * 
  4.  * This file is part of the Ada/Ed-C system.  See the Ada/Ed README file for
  5.  * warranty (none) and distribution info and also the GNU General Public
  6.  * License for more details.
  7.  
  8.  */
  9.  
  10. /* Avoid problems with library routines on PC that require pointer
  11.  * to double instead of double (this is violation of ANSI specs).
  12.  * for now do not compile these modules
  13.  */
  14.  
  15. /* 
  16.  * 11-oct-85    shields
  17.  * fix so can handle most negative integer properly. The changes
  18.  * are more in the nature of a patch than a solution. The proper
  19.  * way is NOT to convert negative values to positive form before
  20.  * processing; instead, positive values should be converted to
  21.  * negative and conversions done on negative values. However, this
  22.  * can be put off until later.
  23.  *
  24.  * 5-jun-85    shields
  25.  * add rat_tol and int_tol which are like rat_toi and int_toi, resp.,
  26.  * except that they return long results. These are needed to support
  27.  * fixed-point in generator. 
  28.  *
  29.  * 1-aug-85    shields
  30.  * add calls to int_free() to free intermediate values known to be dead.
  31.  * add some copies where needed in building new rational values.
  32.  */
  33.  
  34. #include <stdlib.h>
  35. #include <stdio.h>
  36. #include <ctype.h>
  37. #include <string.h>
  38. #include <math.h>
  39. #include "config.h"
  40. #include "miscp.h"
  41.  
  42. typedef struct Rational_s {
  43.     int       *rnum;        /* numerator */
  44.     int       *rden;        /* denominator */
  45. }                Rational_s;
  46. typedef Rational_s * Rational;
  47.  
  48. #include "arithp.h"
  49.  
  50. static int *int_frl(long);
  51. static int *int_ptn(int);
  52. static int *subu_int(int *, int *);
  53. static void int_div(int *, int *, int **, int **);
  54. static int *int__addu(int *, int *);
  55. static int *int__norm(int *);
  56. static int *int_new(int);
  57. static void int_free(int *);
  58. static void rat_free(Rational);
  59. static double pow2(int);
  60.  
  61. /* Some of the procedures want to signal overflow by returning the
  62.  * string 'OVERFLOW'. In C we do this by setting global arith_overflow
  63.  * to non-zero if overflow occurs, zero otherwise
  64.  */
  65. int arith_overflow = 0;
  66.  
  67.  
  68. int    ADA_MIN_INTEGER;
  69. int    ADA_MAX_INTEGER;
  70. int    *ADA_MAX_INTEGER_MP;
  71. int    *ADA_MIN_INTEGER_MP;
  72. long    ADA_MIN_FIXED, ADA_MAX_FIXED;
  73. int    *ADA_MIN_FIXED_MP, *ADA_MAX_FIXED_MP;
  74. /* _LONG form is form that can be held in C (long) */
  75. #ifdef MAX_INTEGER_LONG
  76. long    MIN_INTEGER_LONG;
  77. int    *MAX_INTEGER_LONG_MP;
  78. int    *MIN_INTEGER_LONG_MP;
  79. #endif
  80.  
  81. #define ABS(x) ((x)<0 ? -(x) : (x))
  82. #define SIGN(x) ((x)<0 ? -1 : (x) == 0 ? 0 : 1)
  83.  
  84.  
  85. /*
  86.  * M U L T I P L E     P R E C I S I O N     I N T E G E R
  87.  *
  88.  *       A R I T H M E T I C       P A C K A G E
  89.  *
  90.  *             Robert B. K. Dewar
  91.  *              June 16th, 1980
  92.  *
  93.  * This package of routines implements multiple precision integer
  94.  * arithmetic using what are called the "classical algorithms" in
  95.  * Knuth "Art of Programming", Volume 2, Section 4.3.1. The style
  96.  * of the algorithms follows Knuth fairly closely, and section
  97.  * 4.3.1 can be consulted for further theoretical details.
  98.  *
  99.  * Multiple precision integers are represented as tuples of small
  100.  * integers in the range 0 to BAS - 1, where BAS is a power of 10
  101.  *(actually 10 ** DIGS) which is the base used to represent the
  102.  * long integers. Essentially the successive elements of the tuple
  103.  * are the digits of the representation in base BAS. All integers
  104.  * are normalized so that the first digit is non-zero, except in
  105.  * the case of zero itself. The sign is carried with the first
  106.  * digit value, all remaining digits are always positive.
  107.  *
  108.  * Some examples of the representation, assuming DIGS = 4 and
  109.  * BAS = 10000(note that the choice of BAS as a power of 10
  110.  * is for convenience of the conversion routines, and is not
  111.  * required for correct operation of the arithmetic algorithms).
  112.  
  113.  *     -123456789     [-1 , 2345 , 6789]
  114.  *        0     [0]
  115.  *    123456789     [1 , 2345, 6789]
  116.  
  117.  * The constants BAS and DIGS must be defined as global constants
  118.  * in a program using these routines. It is assumed that the value
  119.  * BAS * BAS - 1 can be represented as a SETL integer value.
  120.  
  121.  * The following routines are provided:
  122.  
  123.  *    int_abs        integer absolute value
  124.  *    int_add        integer addition
  125.  *    int_div        integer division
  126.  *    int_eql        integer equal to
  127.  *    int_exp        raise multiple integer to multiple integer power
  128.  *    int_fri        convert multiple precision integer from integer
  129.  *    int_frs        convert multiple precision integer from string
  130.  *    int_geq        integer greater than or equal to
  131.  *    int_gtr        integer greater than
  132.  *    int_len        number of digits in multiple precision integer
  133.  *    int_leq        integer less than or equal to
  134.  *    int_lss        integer less than
  135.  *    int_mod        integer modulus
  136.  *    int_mul        integer multiplication
  137.  *    int_neq        integer not equal to
  138.  *    int_ptn        integer power of ten
  139.  *    int_quo        integer quotient
  140.  *    int_rem        integer remainder
  141.  *    int_sub        integer subtraction
  142.  *    int_toi        convert integer to C integer
  143.  *    int_tos        integer to string
  144.  *    int_umin    integer unary minus
  145.  */
  146.  
  147. /* Internal procedures have names starting with int__ */
  148.  
  149. int    *int_abs(int *u)                                                /*;int_abs*/
  150. {
  151.     /* Absolute value of multiple precision integer */
  152.     int *pu;
  153.  
  154.     pu = int_copy(u);
  155.     if (pu[1] < 0)
  156.         pu[1] = -pu[1];
  157.     return pu;
  158. }
  159.  
  160. int    *int_add(int *u, int *v)                                        /*;int_add*/
  161. {
  162.     /* Add signed integers
  163.      * This procedure implements the familiar algorithm of comparing
  164.      * the signs, if the signs are the same, then the result is the
  165.      * sum of the magnitudes with the sign of the operands. If the
  166.      * signs differ, then the number with the smaller magnitude is
  167.      * subtracted from the larger magnitude and the result sign is
  168.      * that of the operand with the larger magnitude.
  169.      */
  170.  
  171.     int       *t;
  172.  
  173.     if (u[1] >= 0 && v[1] >= 0)
  174.         return (int__addu(u, v));
  175.     else
  176.         if (u[1] < 0 && v[1] < 0) {
  177.             u[1] = -u[1];
  178.             v[1] = -v[1];
  179.             t = int__addu(u, v);
  180.             u[1] = -u[1];
  181.             v[1] = -v[1];
  182.             t[1] = -t[1];
  183.             return t;
  184.         }
  185.         else {
  186.             int        us, vs;
  187.             if (us = (u[1] < 0)) {
  188.                 u[1] = -u[1];
  189.             }
  190.             if (vs = (v[1] < 0)) {
  191.                 v[1] = -v[1];
  192.             }
  193.             if (int_gtr(u, v)) {
  194.                 t = subu_int(u, v);
  195.                 if (vs) v[1] = -v[1];
  196.                 if (us) {
  197.                     u[1] = -u[1];
  198.                     t[1] = -t[1];
  199.                 }
  200.                 return t;
  201.             }
  202.             else {
  203.                 t = subu_int(v, u);
  204.                 if (us) u[1] = -u[1];
  205.                 if (vs) {
  206.                     v[1] = -v[1];
  207.                     t[1] = -t[1];
  208.                 }
  209.                 return t;
  210.             }
  211.         }
  212. }
  213.  
  214. int    int_eql(int *u, int *v)                                        /*;int_eql*/
  215. {
  216.     /* Compare multiple precision integers for equality */
  217.  
  218.     int        n;
  219.  
  220.     if (u[0] != v[0])
  221.         return FALSE;
  222.  
  223.     for (n = u[0]; n > 0; n--)
  224.         if (u[n] != v[n]) return FALSE;
  225.     return TRUE;
  226. }
  227.  
  228. int    *int_exp(int *u, int *v)                                        /*;int_exp*/
  229. {
  230.     /* Raise a multiple precision integer to a multiple precision integer
  231.      * power using a modified version of the Russian Peasant algorithm
  232.      * for exponentiation.
  233.      */
  234.  
  235.     int        sn;
  236.     int        u1;
  237.     int        i;
  238.     int        carry;
  239.     int        half;
  240.     int       *running, *runningp;
  241.     int       *result, *resultp;
  242.  
  243.     /* Compute sign of result: positive if v is even, the sign of u if
  244.      * v is odd.
  245.      */
  246.  
  247.     sn =((v[v[0]] % 2) == 1) ? SIGN(u[1]) : 1;
  248.     u1 = u[1];
  249.     if (u[1] < 0)
  250.         u[1] = -u1;
  251.     v = int_copy(v);
  252.  
  253.     /* Starting the result at 1 and running at u, loop through the binary
  254.      * digits of v, squaring running each time, and multiplying the result
  255.      * by the current value of running each time a 1-bit is found.
  256.      */
  257.  
  258.     result = int_con(1);
  259.     running = int_copy(u);
  260.  
  261.     while(int_nez(v)) {
  262.         /* If v is odd then multiply result by running. */
  263.  
  264.         if (v[v[0]] % 2 == 1) {
  265.             resultp = result; /* save current value */
  266.             result = int_mul(result, running);
  267.             int_free(resultp); /* free dead value */
  268.         }
  269.  
  270.         /* Square running. */
  271.  
  272.         runningp = running; /* save current value */
  273.         running = int_mul(running, running);
  274.         int_free(runningp); /* free dead value */
  275.  
  276.         /* Halve v. */
  277.  
  278.         carry = 0;
  279.         for (i = 1; i <= v[0]; i++) {
  280.             half = BAS * carry + v[i];
  281.             carry = half % 2;
  282.             v[i] = half / 2;
  283.         }
  284.         v = int__norm(v);
  285.     }
  286.  
  287.     int_free(running);
  288.     int_free(v);
  289.  
  290.     if (sn < 1)
  291.         result[1] = -result[1];
  292.     u[1] = u1;
  293.     return result;
  294. }
  295.  
  296. int    *int_fri(int i)                                                /*;int_fri*/
  297. {
  298.     /* Convert an integer to a multiple precision integer.
  299.      *
  300.      * First check the sign of i.
  301.      */
  302.  
  303.     int        sn = 0;
  304.     int        n = i;
  305.     int       *t;
  306.     int        ti;
  307.     int        d;
  308.  
  309.     if (i < 0) {
  310.         /* Special test for most negative as code below won't work